%% Load the data - Select the Files: 'TurningRate_WildQuasiSteady\Turning_Rates(33-66)','ConditionallySampledTracks'
close all; clear all; clc;
color_option = 't'; % 'arc' for arc, 't' for time
load('.\TurningRate_WildQuasiSteady\Turning_Rates(33-66).mat')
axs = [-0.44 0.44 -0.15 0.3]; %90 degrees
angles = [-90*pi/180 90*pi/180]; % theta = 0 corresponds to parallel to gradient
sigma = 30;
h_size = (6*sigma)-1;
h_filter = fspecial('gaussian', h_size, sigma);
%---9cp
load('.\ConditionallySampledTracks\9cPthetaAndtracks.mat')
for i =1:length(tracks)
    tracks(i).x = imfilter(tracks(i).x,h_filter,'replicate');
    tracks(i).y = imfilter(tracks(i).y,h_filter,'replicate');
end
'Done smoothing tracks'
%
close all;
Fig2b = figure('color',[1 1 1]);
hold on; box on;
colormap jet
len = 160;
[X_left, Y_left, X_right, Y_right] = AverageArclength(tracks, theta, angles, len);
[~, mean_speed] = TrackSpeed(tracks(find(theta(1,:) < (angles(1) + 0.0873) & theta(1,:) > (angles(1) - 0.0873))));
tt = linspace(0,500,1000);
Tau = 30*mean_tau(end); % put it in terms of frames not seconds
conv = 0.65/1000; % mm/pixel
yt = (mean_speed*conv)*(Tau*log(1+exp(2*tt/Tau))-tt-Tau*log(2));
xt1 = 2*Tau*(mean_speed*conv)*(atan(exp(tt/Tau))-atan(1));
xt2 = 2*Tau*(mean_speed*conv)*(atan(-exp(tt/Tau))-atan(-1));
temp = linspace(-1,1,100);
% Black-dotted starting points
startingWidth = 1;
plot(temp*cos(0.0873), temp*sin(0.0873),'--k','linew',startingWidth)
plot(temp*cos(-0.0873), temp*sin(-0.0873),'--k','linew',startingWidth)
x_temp = [];
y_temp = [];
tr = [];
count_1 = 0;
trackWidth = 0.5;
for i = 1:length(tracks)
    if i == 737; continue; end %skip that one pesky track
    if theta(1,i) < (angles(1) + 0.0873) && theta(1,i) > (angles(1) - 0.0873) || theta(1,i) < (angles(2) + 0.0873) && theta(1,i) > (angles(2) - 0.0873)
        x = tracks(i).x(1:len) - tracks(i).x(1);
        y = tracks(i).y(1:len) - tracks(i).y(1);
        z = zeros(len, 1);
        arc = [0; cumsum(sqrt((x(2:end)-x(1:end-1)).^2 + (y(2:end)-y(1:end-1)).^2))];
        t = [1:len]'*mean_speed*conv/0.5; % 0.5 = channel center
        count_1 = count_1 + 1;
        tr(count_1).x = x; tr(count_1).y = y; tr(count_1).t = t;
        % Tracks Plot
        if strcmp(color_option, 'arc')
            surface([x x], [y y], [z z], [arc arc], 'facecol', 'no', 'edgecolor', 'interp','linew',trackWidth)
        else
            surface([x x], [y y], [z z], [t t], 'facecol', 'no', 'edgecolor', 'interp','linew',trackWidth)
        end
    end
end
% Theoretical black curve
theoryWidth = 0.75;
plot(xt1, yt, 'k','linew',theoryWidth)
plot(xt2, yt, 'k','linew',theoryWidth)

arc_left = [0 cumsum(sqrt((X_left(2:end)-X_left(1:end-1)).^2 + (Y_left(2:end)-Y_left(1:end-1)).^2))];
arc_right = [0 cumsum(sqrt((X_right(2:end)-X_right(1:end-1)).^2 + (Y_right(2:end)-Y_right(1:end-1)).^2))];
% Average track: main and border
borderWidth = 3;
meanWidth = 2;
if strcmp(color_option, 'arc')
    plot([flip(X_left(2:len)) X_right(1:len)], [flip(Y_left(2:len)) Y_right(1:len)],'k','linew',borderWidth)
    surface([X_left(1:len)' X_left(1:len)'], [Y_left(1:len)' Y_left(1:len)'],[z(1:len) z(1:len)], [arc_left(1:len)' arc_left(1:len)'], 'facecol', 'no', 'edgecolor', 'interp', 'linew',meanWidth)
    surface([X_right(1:len)' X_right(1:len)'], [Y_right(1:len)' Y_right(1:len)'],[z(1:len) z(1:len)], [arc_right(1:len)' arc_right(1:len)'],'facecol', 'no', 'edgecolor', 'interp', 'linew',meanWidth)
else
    plot([flip(X_left(2:len)) X_right(1:len)], [flip(Y_left(2:len)) Y_right(1:len)],'k','linew',borderWidth)
    surface([X_left(1:len)' X_left(1:len)'], [Y_left(1:len)' Y_left(1:len)'],[z(1:len) z(1:len)], [t(1:len) t(1:len)], 'facecol', 'no', 'edgecolor', 'interp', 'linew',meanWidth)
    surface([X_right(1:len)' X_right(1:len)'], [Y_right(1:len)' Y_right(1:len)'],[z(1:len) z(1:len)], [t(1:len) t(1:len)],'facecol', 'no', 'edgecolor', 'interp', 'linew',meanWidth)
end
daspect([1 1 1])
axis(axs); xticks([]); yticks([]);
set(gca,'color','none')
alw = 0.25; %Border line width
mypapersize = [2.00 1.15]; %[width height]
set(gca,'LineWidth',alw);
set(gca,'XColor',[0 0 0]);set(gca,'YColor',[0 0 0]);
set(gcf, 'PaperUnits', 'inches', 'papersize', mypapersize, 'PaperPosition', [0 0 mypapersize]);
set(Fig2b,'Units','inches','Position',[0 0 mypapersize]);
ax = gca;outerpos = ax.OuterPosition;ti = ax.TightInset;
left = outerpos(1);bottom = outerpos(2);ax_width = outerpos(3);ax_height = outerpos(4);
ax.Position = [left+0.01 bottom+0.035 ax_width-0.02 ax_height-0.04];
% Save Figure
dir = '.\Figure2';
export_fig(Fig2b,[dir,'\9cP_tracks.png'],'-png','-r3000','-transparent');

% Save Data
dir_excel = '.\Source data\MainText\';
tab = 2;
% Mean Trajectory
Mean_X_Left = X_left(1:len)'; Mean_X_Right = X_right(1:len)';
Mean_Y_Left = Y_left(1:len)'; Mean_Y_Right = Y_right(1:len)';
Mean_Time = t(1:len);
writetable(table(Mean_X_Left, Mean_X_Right, Mean_Y_Left, Mean_Y_Right, Mean_Time),[dir_excel,'Figure2.xlsx'],'Sheet',tab,'Range','A3');

% Analytical prediction
Prediction_X_Left = xt1'; Prediction_X_Right = xt2'; Prediction_Y = yt';
writetable(table(Prediction_X_Left,Prediction_X_Right,Prediction_Y),[dir_excel,'Figure2.xlsx'],'Sheet',tab,'Range','G3');

% Tracks
L = 11; % Starting Value
for i = 1:length(tr)
    col = []; col = [char(xlsColNum2Str(L)),'3']; writetable(struct2table(tr(i)),[dir_excel,'Figure2.xlsx'],'Sheet',tab,'Range',col); L = L+3;
end
%
'done save'
% ---contol
load('.\ConditionallySampledTracks\ControlthetaAndtracks.mat')
for i =1:length(tracks)
    tracks(i).x = imfilter(tracks(i).x,h_filter,'replicate'); tracks(i).y = imfilter(tracks(i).y,h_filter,'replicate');
end
'done control track smoothing'
%
close all;
Fig2a = figure('color',[1 1 1]);
hold on; box on; colormap jet
mean_speed_9cp = mean_speed;
[X_left, Y_left, X_right, Y_right] = AverageArclength(tracks, theta, angles, len);
[~, mean_speed] = TrackSpeed(tracks(find(theta(1,:) < (angles(1) + 0.0873) & theta(1,:) > (angles(1) - 0.0873))));
temp = linspace(-1,1,100);
% Black-dotted starting points
plot(temp*cos(0.0873),temp*sin(0.0873),'--k','linew',startingWidth)
plot(temp*cos(-0.0873),temp*sin(-0.0873),'--k','linew',startingWidth)
count_2 = 0;
len = floor(len * mean_speed_9cp/mean_speed); % set the length so that normalized color axis matches
tr_c = [];
for i = 1:length(tracks)
    if theta(1,i) < (angles(1) + 0.0873) && theta(1,i) > (angles(1) - 0.0873) || theta(1,i) < (angles(2) + 0.0873) && theta(1,i) > (angles(2) - 0.0873)
        x = tracks(i).x(1:len) - tracks(i).x(1);
        y = tracks(i).y(1:len) - tracks(i).y(1);
        z = zeros(len, 1);
        arc = [0; cumsum(sqrt((x(2:end)-x(1:end-1)).^2 + (y(2:end)-y(1:end-1)).^2))];
        t = [1:len]'*mean_speed*conv/0.5;
        count_2 = count_2 + 1;
        tr_c(count_2).x = x; tr_c(count_2).y = y; tr_c(count_2).t = t;
        if strcmp(color_option, 'arc')
            surface([x x], [y y], [z z], [arc arc], 'facecol', 'no', 'edgecolor', 'interp','linew',trackWidth)
        else
            surface([x x], [y y], [z z], [t t], 'facecol', 'no', 'edgecolor', 'interp','linew',trackWidth)
        end
    end
end
arc_left = [0 cumsum(sqrt((X_left(2:end)-X_left(1:end-1)).^2 + (Y_left(2:end)-Y_left(1:end-1)).^2))];
arc_right = [0 cumsum(sqrt((X_right(2:end)-X_right(1:end-1)).^2 + (Y_right(2:end)-Y_right(1:end-1)).^2))];
% Average track: main and border
if strcmp(color_option, 'arc')
    plot([flip(X_left(2:len)) X_right(1:len)], [flip(Y_left(2:len)) Y_right(1:len)],'k','linew',borderWidth)
    surface([X_left(1:len)' X_left(1:len)'], [Y_left(1:len)' Y_left(1:len)'],[z(1:len) z(1:len)], [arc_left(1:len)' arc_left(1:len)'], 'facecol', 'no', 'edgecolor', 'interp', 'linew',meanWidth)
    surface([X_right(1:len)' X_right(1:len)'], [Y_right(1:len)' Y_right(1:len)'],[z(1:len) z(1:len)], [arc_right(1:len)' arc_right(1:len)'],'facecol', 'no', 'edgecolor', 'interp', 'linew',meanWidth)
else
    plot([flip(X_left(2:len)) X_right(1:len)], [flip(Y_left(2:len)) Y_right(1:len)],'k','linew',borderWidth)
    surface([X_left(1:len)' X_left(1:len)'], [Y_left(1:len)' Y_left(1:len)'],[z(1:len) z(1:len)], [t(1:len) t(1:len)], 'facecol', 'no', 'edgecolor', 'interp', 'linew',meanWidth)
    surface([X_right(1:len)' X_right(1:len)'], [Y_right(1:len)' Y_right(1:len)'],[z(1:len) z(1:len)], [t(1:len) t(1:len)],'facecol', 'no', 'edgecolor', 'interp', 'linew',meanWidth)
end
st = 0.17; % Line start
plot([0.4; 0.4], [st; st + 0.1022],'-k') % 100 um scale bar: (100/(6.5/10))/1.5051e+03

daspect([1 1 1])
axis(axs)
xticks([]); yticks([]); set(gca,'color','none')
alw = 0.25; % Border line width
mypapersize = [2.00 1.15]; %[width height]
set(gca,'LineWidth',alw);
set(gca,'XColor',[0 0 0]);set(gca,'YColor',[0 0 0]);
set(gcf, 'PaperUnits', 'inches', 'papersize', mypapersize, 'PaperPosition', [0 0 mypapersize]);
set(Fig2a,'Units','inches','Position',[0 0 mypapersize]);

ax = gca;outerpos = ax.OuterPosition;ti = ax.TightInset;
left = outerpos(1);bottom = outerpos(2);ax_width = outerpos(3);ax_height = outerpos(4);
ax.Position = [left+0.01 bottom+0.035 ax_width-0.02 ax_height-0.04];

% Save Data
clc;
dir_excel = '.\Source data\MainText\';
tab = 1;
% Mean Trajectory
Mean_X_Left = X_left(1:len)'; Mean_X_Right = X_right(1:len)';
Mean_Y_Left = Y_left(1:len)'; Mean_Y_Right = Y_right(1:len)';
Mean_Time = t(1:len);
writetable(table(Mean_X_Left, Mean_X_Right, Mean_Y_Left, Mean_Y_Right, Mean_Time),[dir_excel,'Figure2.xlsx'],'Sheet',tab,'Range','A3');
% Tracks
L = 7; % Starting Value
for i = 1:length(tr_c)
    col = []; col = [char(xlsColNum2Str(L)),'3']; writetable(struct2table(tr_c(i)),[dir_excel,'Figure2.xlsx'],'Sheet',tab,'Range',col); L = L+3;
end
% Save Figure
dir = '.\Figure2'
export_fig(Fig2a,[dir,'\Control_curvetracks_highres.png'],'-png','-r3000','-transparent')
'Done Save: Control'
% Colorbar
close all;
ColorbarFig = figure('color',[1 1 1]); hold on; box on;
colormap jet

% Average track: main and border
borderWidth = 3;
meanWidth = 2;
if strcmp(color_option, 'arc')
    plot([flip(X_left(2:len)) X_right(1:len)], [flip(Y_left(2:len)) Y_right(1:len)],'k','linew',borderWidth)
    surface([X_left(1:len)' X_left(1:len)'], [Y_left(1:len)' Y_left(1:len)'],[z(1:len) z(1:len)], [arc_left(1:len)' arc_left(1:len)'], 'facecol', 'no', 'edgecolor', 'interp', 'linew',meanWidth)
    surface([X_right(1:len)' X_right(1:len)'], [Y_right(1:len)' Y_right(1:len)'],[z(1:len) z(1:len)], [arc_right(1:len)' arc_right(1:len)'],'facecol', 'no', 'edgecolor', 'interp', 'linew',meanWidth)
else
    plot([flip(X_left(2:len)) X_right(1:len)], [flip(Y_left(2:len)) Y_right(1:len)],'k','linew',borderWidth)
    surface([X_left(1:len)' X_left(1:len)'], [Y_left(1:len)' Y_left(1:len)'],[z(1:len) z(1:len)], [t(1:len) t(1:len)], 'facecol', 'no', 'edgecolor', 'interp', 'linew',meanWidth)
    surface([X_right(1:len)' X_right(1:len)'], [Y_right(1:len)' Y_right(1:len)'],[z(1:len) z(1:len)], [t(1:len) t(1:len)],'facecol', 'no', 'edgecolor', 'interp', 'linew',meanWidth)
end
xticks([]); yticks([]);
cbar = colorbar;
axis off
cbar.LineWidth = 0.5; cbar.Color = 'k'; cbar.TickLength = 0.02; cbar.Ticks = [0.1 0.2 0.3 0.4 0.5]; cbar.TickLabels = [];
set(gca,'color','none')

alw = 0.5; % Border line width
mypapersize = [1.8031 1.1331]; %[width height]
 
set(gca,'LineWidth',alw);
set(gca,'XColor',[0 0 0]);set(gca,'YColor',[0 0 0]);
set(gcf, 'PaperUnits', 'inches', 'papersize', mypapersize, 'PaperPosition', [0 0 mypapersize]);
set(ColorbarFig,'Units','inches','Position',[0 0 mypapersize]);

dir = '.\Figure2'
print(ColorbarFig, '-painters','-dpdf', '-r1200', [dir,'\ColorBar.pdf']); %r1200
%% FUNCTION: Getting average trajectories in time from set angles
function [X_left, Y_left, X_right, Y_right] = AverageArclength(tracks, theta, angles, len)
X_left = []; Y_left = [];X_right = []; Y_right = [];
ind_left = find(theta(1,:) <(angles(1) + 0.0873) & theta(1,:) > (angles(1) - 0.0873));
ind_right = find(theta(1,:) <(angles(2) + 0.0873) & theta(1,:) > (angles(2) - 0.0873));
for i = 1:len
    temp_X_left =[]; temp_Y_left = []; temp_X_right = []; temp_Y_right = [];
    for ii = 1:length(tracks)
    	if isnan(tracks(ii).x(i))
    		continue
    	end
        if ~isempty(find(ind_left == ii,1))
            temp_X_left = [temp_X_left, tracks(ii).x(i) - tracks(ii).x(1)];
            temp_Y_left = [temp_Y_left, tracks(ii).y(i) - tracks(ii).y(1)];
        elseif ~isempty(find(ind_right == ii,1))
            temp_X_right = [temp_X_right, tracks(ii).x(i) - tracks(ii).x(1)];
            temp_Y_right = [temp_Y_right, tracks(ii).y(i) - tracks(ii).y(1)];
        end
    end
    X_left(i) = nanmean(temp_X_left); Y_left(i) = nanmean(temp_Y_left);
    X_right(i) = nanmean(temp_X_right); Y_right(i) = nanmean(temp_Y_right);
end
end

%% FUNCTION: speed from tracks
function [speed, mean_speed] = TrackSpeed(tracks)
try
    u = [tracks.u]; v = [tracks.v];
catch
    u = vertcat(tracks.u); v = vertcat(tracks.v);
end
speed = sqrt(u.^2 + v.^2); mean_speed = nanmean(speed);

end
%% FUNCTION: Excel Column Labels
function [colChar]=xlsColNum2Str(colNum)
% Reference: Matt G (2021). Excel Column Number To Column Name (https://www.mathworks.com/matlabcentral/fileexchange/15748-excel-column-number-to-column-name), MATLAB Central File Exchange. Retrieved March 23, 2021. 
%XLSCOLNUM2STR takes in an array of numbers and returns a cellular array
%of the same size with cell of corresponding Excel column names.
%
%For example:
%n=[1  10;
%   53 256]
%c=xlsColNum2Str(n);
%c={'A' , 'J';
%   'BA', 'IV'}
%Note: up to Excel 2003 the number of columns was limited to 256, as of
%Excel 2007 the number of columns has increased to 16,384 or 'XFD'
%This function is designed to take accept any integer so proper handling 
%of the number of columns should be taken care of outside this function
    colChar=cell(size(colNum)); %blank cell array
    
    % find max number of characters (AA n=2)
    numOfChars=ceil(max(colNum)/26)-1;
    n=1;
    
    while numOfChars>=1
        numOfChars=ceil(numOfChars/26)-1;
        n=n+1;
    end    
    
    remainder=num2cell(colNum);
    
    for s=n:-1:1
        if s>1
            %find limits
            % if n=2 then the columns go from AA to ZZ or 27 to 702
            L=sum(26.^(1:s-1))+1; % lower limit
            U=sum(26.^(1:s));   %upper limit
            %place current character to right of previous
            colChar(colNum>=L & colNum<=U) = ...
                cellfun(@(x,y) ([x char(ceil((y-(L-1))/26^(s-1))+64)]),...
                colChar(colNum>=L & colNum<=U),...      % x
                remainder(colNum>=L & colNum<=U),...    % y
                'UniformOutput',false);
            %calculate the remaining string
            %for example if last string was 'ABA' the 'A' was placed to the
            %right of the previous string and now 'BA' is remaining
            remainder(colNum>=L & colNum<=U)=...
                cellfun(@(x,y) (y-26^(s-1)*(double(x(end))-64)),...
                colChar(colNum>=L & colNum<=U),...
                remainder(colNum>=L & colNum<=U),'UniformOutput',false);
            colNum=cell2mat(remainder);
        else
             colChar=cellfun(@(x,y) ([x char(y+64)]),...
                 colChar,remainder,'UniformOutput',false);
        end
    end    
end